The amount of techniques for processing satellite imagery can be overwhelming. This tutorial is an attempt to summarize and organize some of these techniques, with a special focus on those available in R (either native or via third-party software, e. g. SAGA). Other methods not available in R are also described, but just enough to provide a greater context and be aware about the strengths and limitations using R as a tool for pre-processing.
The topic of this tutorial is pansharpening , a group of techniques that increase the spatial and spectral resolution of satellite images by fusing images with complementary features. The tutorial gives: An introduction to the basics of pansharpening (Section 1), an overview of methods (Section 2), and descriptions of some methods (Section 3), including a brief theoretical explanations and simple and reproducible examples. To follow this tutorial, install RStoolbox(Leutner et al., 2019), tmap (Tennekes, 2018) and raster (Hijamns 2020), if you don’t have them already:
install.packages("raster")
install.packages("RStoolbox")
install.packages("tmap")To load them:
library("raster")## Loading required package: splibrary("RStoolbox")## Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
## when loading 'dplyr'library("tmap")For convenience, we built the compare_vis() function for this tutorial. The function has two inputs, the original and the pansharpened multi-spectral images. The function is a wrapper of tmap and will display side-by-side interactive maps of RGB images. The function is available in this repository, and can be loaded into R as follows:
source("./R/compare_vis.R")Pansharpening methods generally fuse a panchromatic and a multi-spectral (or hyper-spectral) image:
A panchromatic (PAN) image is an image with greater spatial details (but little spectral information). Instruments retrieving PAN images sense broader bandwidths (visible and to some degree near-infrared) to increase the signal-to-noise ratio and capture finer spatial details. This also means that panchromatic images have a single band and they do not distinguish intensities from discrete wavelengths. In this tutorial, we use the PAN image under the path ./Data/small_sample, with a resolution of \(15 \times 15m^2\) resolution:
ms.file <- "./Data/small_sample/pan_example.tif"
pan.img <- raster(ms.file)
pan.img## class : RasterLayer
## dimensions : 620, 799, 495380 (nrow, ncol, ncell)
## resolution : 15.00038, 14.99936 (x, y)
## extent : 605264.5, 617249.8, 4736529, 4745828 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : D:/github/pansharpen_r/Data/small_sample/pan_example.tif
## names : pan_example
## values : 0.1180896, 0.9924163 (min, max)Multi(hyper)-spectral (MS) images carry information about the radiance/reflectance from narrower wavelengths. The difference between multi and hyper-spectral lays on the number of bands (tens vs. hundreds) and how narrow the bandwidths are. Our multi-spectral image has 4 bands; red, green, blue, and near-infrared (NIR). The image can be found under ./Data/small_sample and has a resolution of \(30 \times 30m^2\) :
ms.img <- stack("./Data/small_sample/ms_example.tif")
ms.img## class : RasterStack
## dimensions : 310, 400, 124000, 4 (nrow, ncol, ncell, nlayers)
## resolution : 29.9976, 30.00156 (x, y)
## extent : 605257.9, 617256.9, 4736523, 4745823 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : ms_example.1, ms_example.2, ms_example.3, ms_example.4
## min values : 0.1105364, 0.1246357, 0.1405966, 0.1008164
## max values : 0.7430533, 0.6397192, 0.4251927, 0.9167315PAN and MS images were captured by Landsat-8, on June 30st 2015. They both belong to the Collection 1, and processing Level 1. Before pansharpening, the literature recommends the co-registration of the PAN and MS so they are properly aligned. This is true even when if they belong to the same platform (Wegmann et al., 2016). The function coregisterImages() from RStoolbox applies an automated co-registration algorithm. It shifts a slave image along the vertical and horizontal axes over a reference image to maximize the mutual information. This function demands both master and slave to have the same number of bands. It wasn’t conceived to co-register MS and PAN images. Here we run an informal test with a single band from MS and the PAN, to illustrate the application of this function:
# shift from -10 to +10 pixels horizontally and vertically
img.shift <- coregisterImages(ms.img[[1]],pan.img,shift=10,reportStats=TRUE)
img.shift$bestShift## x y
## 221 0 0As the best shift is \(0,0\), we interprete that PAN and MS images are well aligned along both axes.
Pansharpening methods are classified in two main categories (Alparone et al., 2015):
Component Substitution (CS) methods: These methods project the multi-spectral image into a new vector space. The new vector space disentangles the spatial structure from the spectral information. The component representing the spatial structure is replaced by the PAN image. Finally, the image is projected back to the original vector space. Among the most popular CS techniques are the Intensity-Hue-Saturation method (IHS), the Principal Component Analysis (PCA), and the Gram-Schmidt (GS) spectral sharpening.
Multi-Resolution Analyses (MRA): This group of methods extract the spatial structure from the PAN image through multi-resolution decomposition. The spatial details are then injected into an interpolated version of the multi-spectral image. The most popular are the (decimated) Discrete Wavelet Transform (DWT), Undecimated Wavelet Transform (UDWT), A-Trous Wavelet Transform (ATW), and the Laplacian Pyramid (LP).
A third group combines both CS and MRA, and thus they are called hybrid. Current software in (or linked to) R covers some CS methods. They provide a high geometrical quality and they are simple and computationally fast. However, they do not consider potential discrepancies between the PAN and the MS image, leading, in some circumstances, to spectral distortions. They should be used with care (or avoided) for the analysis of spectral signatures (Alparone et al., 2015; Wegmann et al., 2016).
Mains steps in CS sharpening are: (1) interpolate the MS image to match the scale of the PAN image, (2) define and apply spectral weights to each band to obtain the intensity image (i.e., an image representing the darkness-lightness of colors), (3) match the histograms of the PAN and the intensity image, so they show similar frequencies for each intensity value, (4) calculate and sum the injection gains to the original image (i.e. the proportional part that goes to each band of the original image). Variations of CS differ in the definition of the spectral weights in step 2 and the calculation of the injection gains in step 4.
The RStoolbox package implements three variants of the CS (Leutner et al., 2019) with the panSharpen() function; the Brovey Transform (BT), the Intensity-Hue-Saturation (IHS), and the Principal Component Analysis (PCA). Another alternative is using the image processing module from SAGA , which additionally includes the Spectral Sharpening. We will check how SAGA may boost the computational performance of pansharpening when dealing with large images. Due to the need of third-party software and the methodological overlap between RStoolbox and SAGA, we relegate its use to the end of this tutorial.
BT (Gillespie et al., 1987) is one of the simplest pansharpening techniques. BT applies ratios between the bands of the MS and the PAN to obtain the sharpened image. More specifically, panSharpen() first interpolates the MS image (\(\tilde{MS}\)) to match the resolution of the PAN image using nearest neighbors, then calculates the intensity image as the mean of all bands (\(I = \frac{1}{N} \sum_{k = 1}^{N}{MS_k}\)), and finally injects the details to the original image by multiplying ratio between the PAN and intensity image:
\(MS_{PAN,k}=\tilde{MS,k}\dot\frac{PAN}{I}\)
The Brovey transform can be applied setting the argument method = "brovey". BT only deals with the Red-Green-Blue (RGB) bands, so we must specify which layer in the MS image corresponds to R, G, and B. The red, green, and blue bands in ms.img are the 1st, 2nd, and 3rd layers respectively. Thus, we run:
library("RStoolbox")
bro.img <- panSharpen(ms.img, pan.img,r=1,g=2,b=3,method="brovey")We can visually compare the original and the sharpened image with the function compare_vis() :
compare_vis(bro.img, ms.img, titles = c("Brovey", "Original"))## tmap mode set to interactive viewing## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1IHS (Carper et al., 1990) begins with transformation of the interpolated RGB image (\(\tilde{MS}\)) into the Intensity-Hue-Saturation color space, which is a model of visual perception of colors. With this decomposition we extract the intensity image as a linear combination of the RGB bands (generally expressed as \(I=\sum_{i=0}^{N}w_i\tilde{MS_i}\)). The histograms from the intensity and the PAN image are matched. Then, the PAN replaces the intensity band, and the IHS images is transformed back to Red-Green-Blue. In addition to
ihs.img <- panSharpen(ms.img, pan.img,r=1,g=2,b=3,method="ihs")Let’s display the original and the pansharpened image:
compare_vis(ihs.img, ms.img = ms.img, titles = c("IHS", "Original"))## tmap mode set to interactive viewingNote that the classical approach only deals with the Red-Green-Blue spectral bands (Wegmann et al., 2016). It preserves the color balance relatively well (Wegmann et al., 2016).
PCA (Licciardi et al., 2012) transforms the spectral bands such that they are projected onto new axes, called Principal Components, which are statistically uncorrelated with each other. The components are generally sorted with decreasing variance and the frequent assumption is that the 1st component (PC1) is the one that spectrally matches the PAN. The weights to compute the intensity image (step 2) and injection gains (step 4) are derived from the PCA.
pca.img <- panSharpen(ms.img, pan.img,method="pca")Let’s have a look at the results:
compare_vis(pca.img, ms.img = ms.img, titles = c("PCA", "Original"))## tmap mode set to interactive viewingThe strengths of the PCA are the high fidelity reproducing sharp edges and and the ability to deal with more than three bands. The weakness is that it can lead to severe color distortions (Wegmann et al., 2016), as shown in this case. The success of the fusion will depend on which extent the assumption \(PC1 \sim PAN\) is fulfilled (Alparone et al., 2015).
As a final result, we display all the images together:
compare_vis(bro.img, ihs.img, pca.img, ms.img, titles = c("BT", "IHS", "PCA", "Original"), lyout = c(1,4))## tmap mode set to interactive viewingThe multi-spectral images from above are purposely small (\(310 \times400\) pixels). R can be slow when dealing with large images. To illustrate the issue, we repeat the analysis with larger PAN and MS images (\(1579 \times 1882\) pixels for the MS image). Both images are saved in the ./Data/large_sample directory. Load them in R as follows:
ms.lrg <- stack("./Data/large_sample/ms_large_example.tif")
pn.lrg <- raster("./Data/large_sample/pan_large_example.tif")The function mark() from the bench package (Hester, 2020) tracks the running times in R . The argument time_unit = "m" ensures consistent time units (here minutes) across tests. We run the three pansharpening techniques from above using the large images. Please, be patient as next lines may take a few minutes)
library(bench)
bro.rstl <- mark(
panSharpen(ms.lrg, pn.lrg,r=1,g=2,b=3,method="brovey"),
time_unit = "m")## Warning: Some expressions had a GC in every iteration; so filtering is disabled.bro.rstl$min## [1] 1.883032ihs.rstl <- mark(
panSharpen(ms.lrg, pn.lrg,r=1,g=2,b=3,method="ihs"),
time_unit = "m")## Warning: Some expressions had a GC in every iteration; so filtering is disabled.ihs.rstl$min## [1] 1.327017pca.rstl <- mark(
panSharpen(ms.lrg, pn.lrg, method="pca"),
time_unit = "m")## Warning: Some expressions had a GC in every iteration; so filtering is disabled.pca.rstl$min## [1] 2.011058The code takes \(1.89\), \(1.37\), and \(1.85\) minutes to pansharpen the large MS image using an Intel(R) Core(TM) i7-4710HQ CPU at 2.50 GHz. Under similar (or worse) circumstances, it might be worthwhile to consider external libraries such as SAGA. This GIS software provides pansharpening techniques within the module called imagery tools. SAGA offers compiled code, which is significantly faster than R. Before continuing with this tutorial, make sure SAGA is installed. Once installed, find the saga_cmd.exe file and paste the path in the following command. For me:
saga.path <- normalizePath("C:\\OSGeo4W64\\apps\\saga-ltr\\saga_cmd.exe")The github repository provides custom functions to interact with the pansharpening libraries in SAGA. The functions are inside the ./R folder and there is one function for each pansharpening technique; the Brovey, IHS, and PCA methods:
source("./R/pansharpen_brovey.R")
source("./R/pansharpen_ihs.R")
source("./R/pansharpen_pca.R")Then, the same analysis can be performed with SAGA as follows:
bro.saga <- mark(pansharpen_brovey(saga.path, ms.lrg, pn.lrg), time_unit = "m")## Warning in dir.create(tmp.dir): 'D:\github\pansharpen_r\tmp' already existsihs.saga <- mark(pansharpen_ihs(saga.path, ms.lrg, pn.lrg), time_unit = "m")
pca.saga <- mark(pansharpen_pca(saga.path, ms.lrg, pn.lrg), time_unit = "m")Each of the pansharpening techniques takes \(0.48\), \(0.50\), and \(0.7\) minutes, which saves around \(25-50\%\) of the time relative to RStoolbox. Note that there is still room for improvement in the custom functions to achieve better computational efficiency. You may also consider using specialized packages, such as RSAGA(Brenning, 2018), to achieve better communications between R and SAGA. Also, using SAGA from R does not always means a faster implementation. With small images, the cost of calling SAGA is greater than the savings from compiled code and native R might be preferred.
Alparone, L., Aiazzi, B., Baronti, S., & Garzelli, A. (2015). Pansharpening of Multispectral Images, in Remote Sensing Image Fusion. Crc Press.
Brenning A., Bangs D., and Becker M. (2018). RSAGA: SAGA Geoprocessing and Terrain Analysis. R package version 1.3.0. (URL: https://CRAN.R-project.org/package=RSAGA)
Carper W., Lillesand T., and Kiefer , R. (1990). The use of intensity-hue-saturation transformations for merging SPOT panchromatic and multispectral image data. Photogrammetric Engineering and Remote Sensing, 56(4):459–467.
Gillespie, A. R., Kahle, A. B., & Walker, R. E. (1987). Color enhancement of highly correlated images. II. Channel ratio and “chromaticity” transformation techniques. Remote Sensing of Environment, 22(3), 343-365. (URL: https://doi.org/10.1016/0034-4257(87)90088-5)
Hester J. (2020). bench: High Precision Timing of R Expressions. R package version 1.1.1. (URL: https://CRAN.R-project.org/package=bench)
Hijmans R. (2020). raster: Geographic Data Analysis and Modeling. R package version 3.3-13. (URL: https://CRAN.R-project.org/package=raster)
Leutner B., Horning N., & Schwalb-Willmann, J. (2019). RStoolbox: Tools for Remote Sensing Data Analysis. R package version 0.2.6. (URL: https://CRAN.R-project.org/package=RStoolbox)
Licciardi, G. A., Khan, M. M., Chanussot, J., Montanvert, A., Condat, L., & Jutten, C. (2012). Fusion of hyperspectral and panchromatic images using multiresolution analysis and nonlinear PCA band reduction. EURASIP Journal on Advances in Signal processing, 2012(1), 1-17. (URL: [https://doi.org/](https://doi.org/10.1016/0034-4257(87)90088-5)10.1109/IGARSS.2011.6049466)
Tennekes, M. (2018). tmap: Thematic Maps in R. Journal of Statistical Software 84 (6), 1-39. (URL: https://doi.org/10.18637/jss.v084.i06)..)
Wegmann, M., Leutner, B., & Dech, S. (Eds.). (2016). Remote sensing and GIS for ecologists: using open source software. Pelagic Publishing Ltd.